function Lrk4_error

%  Plots the error at t=tmax using the C-N & L-RK4 methods
%    as a function of the number of time points used, and
%  for the heat equation problem
%       diff(u,x,x) = diff(u,t)   for xL < x < xr, 0 < t < tmax
%  where
%      u = 0  at x=xL,xR  and  u = sin(2*pi*x) at t = 0

% clear all previous variables and plots
clear *
clf

% set parameters
tmax=0.1;
xL=0;
xR=1;

N=20;

% generate the points along the x-axis, x(1)=xL and x(N+2)=xR
x=linspace(xL,xR,N+2);
h=x(2)-x(1);

ii=0;  ff=1;  MM=2;
for i=1:21
	if i==16
		ff=5;
	elseif i==21
		ff=10;
	elseif i==32
		ff=100;
	end;
	MM=MM+ff;
	ii=ii+1; points1(ii)=MM-1;
	t=linspace(0,tmax,MM);
	k=t(2)-t(1);
	lamda=k/h^2;
	errorC1(ii)=errorC(x,t,N+2,MM,tmax,lamda);
end;	

ii=0;  ff=1;  MM=54;
for i=1:12
	if i==5
		ff=10;
	elseif i==9
		ff=20;
	end;
	MM=MM+ff;
	ii=ii+1; points2(ii)=MM-1;
	t=linspace(0,tmax,MM);
	k=t(2)-t(1);
	lamda=k/h^2;
	errorRK1(ii)=errorRK(x,t,N+2,MM,tmax,lamda);
end;	


N=40;

% generate the points along the x-axis, x(1)=xL and x(N+2)=xR
x=linspace(xL,xR,N+2);
h=x(2)-x(1);

ii=0; ff=2; MM=3;
for i=1:25
	if i==10
		ff=1;
	elseif i==16
		ff=3;
	elseif i==20
		ff=15;
	elseif i==24
		ff=50;
	end;
	MM=MM+ff;
	ii=ii+1; points3(ii)=MM-1;
	t=linspace(0,tmax,MM);
	k=t(2)-t(1);
	lamda=k/h^2;
	errorC2(ii)=errorC(x,t,N+2,MM,tmax,lamda);
end;	

ii=0; ff=1; MM=232;
% min MM=233
for i=1:11
	if i==5
		ff=100;
	end;
	MM=MM+ff;
	ii=ii+1; points4(ii)=MM-1;
	t=linspace(0,tmax,MM);
	k=t(2)-t(1);
	lamda=k/h^2;
	errorRK2(ii)=errorRK(x,t,N+2,MM,tmax,lamda);
end;	

% get(gcf)
%set(gcf,'Position', [1260 559 595 335]);
plotsize(595, 335)

% plot results
loglog(points1,errorC1,'--bo','MarkerSize',7)
hold on
loglog(points2,errorRK1,'-bd','MarkerSize',7)
loglog(points3,errorC2,'--rs','MarkerSize',7)
loglog(points4,errorRK2,'-r*','MarkerSize',7)
legend(' N = 20 (CN)',' N = 20 (RK4)',' N = 40 (CN)',' N = 40 (RK4)',3);

% axes used in plot
xlabel('Number of Time Points','FontSize',14,'FontWeight','bold')
ylabel('Error','FontSize',14,'FontWeight','bold')

% have MATLAB use certain plot options (all are optional)
grid on
set(gca,'MinorGridLineStyle','none')
set(gca,'FontSize',14)
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');
hold off
say=['Comparison of errors when u(x,0)=sin(2\pix)'];
title(say,'FontSize',14,'FontWeight','bold')


% error calculation for c-n
function q=errorC(x,t,N,M,tmax,lamda)
for i=1:N
	exact(i)=exp(-4*pi*pi*tmax)*sin(2*pi*x(i));
end;
h=x(2)-x(1);  k=t(2)-t(1);
ue=crank(x,t,N,M,h,k,lamda);
q=norm(ue(:,M)-exact',inf);


% error calculation for rk4 '
function q=errorRK(x,t,N,M,tmax,lamda)
for i=1:N
	exact(i)=exp(-4*pi*pi*tmax)*sin(2*pi*x(i));
end;
h=x(2)-x(1);  k=t(2)-t(1);
ue=rk4(x,t,N,M,h,k);
q=norm(ue(:,M)-exact',inf);


% c-n method '
function UC=crank(x,t,N,M,h,k,lamda)
UC=zeros(N,M);
for i=1:N
	UC(i,1)=g(x(i));
end;
a=-2*(1+lamda)*ones(1,N); b=lamda*ones(1,N); c=b; z=zeros(1,N);
a(1)=1; c(1)=0; a(N)=1; b(N)=0;
for j=2:M
	for i=2:N-1
		z(i)=-lamda*UC(i+1,j-1)-2*(1-lamda)*UC(i,j-1)-lamda*UC(i-1,j-1)+k*f(x(i),t(j))+k*f(x(i+1),t(j-1));
	end;
	UC(:,j) = tridiag( a, b, c, z )';
end;

% rk4 method '
function UC=rk4(x,t,N,M,h,k)

UC=zeros(N,M);
for i=1:N
	UC(i,1)=g(x(i));
end;
N2=N-2;
alpha=1/(h*h);
A=diag(-2*alpha*ones(N2,1))+diag(alpha*ones(N2-1,1),1)+diag(alpha*ones(N2-1,1),-1);
y=zeros(N2,1);
for i=1:N2
	y(i)=g(x(i+1));
end;
for j=2:M
	k1=k*A*y;
	k2=k*A*(y+0.5*k1);
	k3=k*A*(y+0.5*k2);
	k4=k*A*(y+k3);
	y=y+(k1+2*k2+2*k3+k4)/6;
	UC(:,j) = [0; y; 0];
end;


% subfunction f(x,t)
function q=f(x,t)
q=0;

% subfunction g(x)
function q=g(x)
q=sin(2*pi*x);

% tridiagonal solver
function y = tridiag( a, b, c, f )
N = length(f);
v = zeros(1,N);   
y = v;
w = a(1);
y(1) = f(1)/w;
for i=2:N
    v(i-1) = c(i-1)/w;
    w = a(i) - b(i)*v(i-1);
    y(i) = ( f(i) - b(i)*y(i-1) )/w;
end;
for j=N-1:-1:1
   y(j) = y(j) - v(j)*y(j+1);
end;

% subfunction plotsize
function plotsize(width,height)
siz=get(0,'ScreenSize');
bottom=max(siz(4)-height-95,1);
set(gcf,'Position', [2 bottom width height]);